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A Molecular Dynamics (MD) parallel to the Control Volume (CV) formulation of fluid mechanics is devel- 
oped by integrating the formulas of Irving and Kirkwood, J. Chem. Phys. 18, 817 (1950) over a finite cubic 
volume of molecular dimensions. The Lagrangian molecular system is expressed in terms of an Eulerian CV, 
which yields an equivalent to Reynolds' Transport Theorem for the discrete system. This approach casts the dy- 
namics of the molecular system into a form that can be readily compared to the continuum equations. The MD 
equations of motion are reinterpreted in terms of a Lagrangian-to-Control- Volume {CCV) conversion function 
$i, for each molecule i. The CCV function and its spatial derivatives are used to express fluxes and relevant 
forces across the control surfaces. The relationship between the local pressures computed using the Volume 
Average (VA, Lutsko, J. Appl. Phys 64, 1 152 (1988) ) techniques and the Method of Planes (MOP , Todd et al, 
Phys. Rev. E 52, 1627 (1995) ) emerges naturally from the treatment. Numerical experiments using the MD CV 
method are reported for equilibrium and non-equilibrium (start-up Couette flow) model liquids, which demon- 
strate the advantages of the formulation. The CV formulation of the MD is shown to be exactly conservative, 
and is therefore ideally suited to obtain macroscopic properties from a discrete system. 
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I. INTRODUCTION 

The macroscopic and microscopic descriptions of mechan- 
ics have traditionally been studied independently. The former 
invokes a continuum assumption, and aims to reproduce the 
large-scale behaviour of solids and fluids, without the need to 
resolve the micro-scale details. On the other hand, molecu- 
lar simulation predicts the evolution of individual, but inter- 
acting, molecules, which has application in nano and micro- 
scale systems. Bridging these scales requires a mesoscopic 
description, which represents the evolution of the average of 
many microscopic trajectories through phase space. It is ad- 
vantageous to cast the fluid dynamics equations in a consis- 
tent form for both the molecular, mesoscale and continuum 
approaches. The current works seeks to achieve this objec- 
tive by introducing a Control Volume (CV) formulation for 
the molecular system. 

The Control Volume approach is widely adopted in con- 
tinuum fluid mechanics, where Reynolds Transport Theorem 
|Q]] relates Newton's laws of motion for macroscopic fluid 
parcels to fluxes through a CV. In this form, fluid mechanics 
has had great success in simulating both fundamental yj, |3|] 
and practical 10-Sl flows. However, when the continuum as- 
sumption fails, or when macroscopic constitutive equations 
are lacking, a molecular-scale description is required. Exam- 
ples include nano-flows, moving contact lines, solid-liquid 
boundaries, non-equilibrium fluids, and evaluation of trans- 
port properties such as viscosity and heat conductivity H . 

Molecular Dynamics (MD) involves solving Newton's 
equations of motion for an assembly of interacting discrete 
molecules. Averaging is required in order to compute proper- 
ties of interest, e.g. temperature, density, pressure and stress, 
which can vary on a local scale especially out of equilib- 
rium Q. A rigorous link between mesoscopic and continuum 
properties was established in the seminal work of Irving and 
Kirkwood J8Q, who related the mesoscopic Liouville equa- 



PACS number(s): 05.20.y, 47.11.Mn, 31.15.xv 



tion to the differential form of continuum fluid mechanics. 
However, the resulting equations at a point were expressed 
in terms of the Dirac 8 function — a form which is difficult 
to manipulate and cannot be applied directly in a molecular 
simulation. Furthermore, a Taylor series expansion of the 
Dirac S functions was required to express the pressure ten- 
sor. The final expression forpressure tensor is neither easy 
to interpret nor to compute H. As a result, there have been 
numerous attempts to develop an ex pres sion for the pressure 
tensor for use in MD simulation ll9 H2lll . Some of these ex- 
pressions have been shown to be equivalent in the appropriate 
limit. For example, Heyes et al. 12211 ) demonstrated equiva- 
lence between Method of Planes (MOP Todd et al. (H]) and 
Volume Average (VA Lutsko I16IP at a surface. 

In order to avoid use of the Dirac 5 function, the current 
work adopts a Control Volume representation of the MD sys- 
tem, written in terms of fluxes and surface stresses. This ap- 
proach is in part motivated by the success of the control vol- 
ume formulation in continuum fluid mechanics. At a molecu- 
lar scale, control volume analyses of NEMD simulations can 
facilitate evaluation of local fluid properties. Furthermore, 
the C V method also lends itself to coupling schemes between 
the continuum and molecular descriptions IE3H34I1 . 

The equations of continuum fluid mechanics are presented 
in Section ITLAl followed by a review of the Irving and Kirk- 
wood [80 procedure for linking continuum and mesoscopic 
properties in Section HTBl In sectionUm a Lagrangian to Con- 
trol Volume {CCV) conversion function is used to express the 
mesoscopic equations for mass and momentum fluxes. Sec- 
tion IIIICI focuses on the stress tensor, and relates the cur- 
rent formulation to established definitions within the litera- 
ture llo[m[l7ll . In Section lTVl the CV equations are derived 
for a single microscopic system, and subsequently integrated 
in time in order to obtain a form which can be applied in MD 
simulations. The conservation properties of the CV formula- 
tion are demonstrated in NEMD simulations of Couette flow 
in Section llVCl 
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II. BACKGROUND 

This section summarizes the theoretical background. First, 
the macroscopic continuum equations are introduced, fol- 
lowed by the mesoscopic equations which describe the evolu- 
tion of an ensemble average of systems of discrete molecules. 
The link between the two descriptions is subsequently dis- 
cussed. 

A. Macroscopic Continuum Equations 

The continuum conservation of mass and momentum bal- 
ance can be derived in an Eulerian frame by considering the 
fluxes through a Control Volume (CV). The mass continuity 
equation can be expressed as, 



!£*" 



pu ■ dS, 



(1) 



where p is the mass density and u is the fluid velocity. The 
rate of change of momentum is determined by the balance of 
forces on the CV, 



a_ 

m 



pudV - 



puu ■ dS + F surface + F body . (2) 



The forces are split into ones which act on the bounding sur- 
faces, F sur f ace , and body forces, Fbody Surface forces are 
expressed in terms the pressure tensor, II, on the CV sur- 
faces, 



and rj, respectively, and N is the number of molecules in a 
single system. The momentum density at a point in space is 
similarly defined by, 



A' 



p(r, t)u(r, t) = J2{ Pi S ( r i - r );/ 



t=l 



(8) 



where the molecular momentum, p, ; = mjTj. Note that Pj is 
the momentum in the laboratory frame, and not the peculiar 
value Pj which excludes the macroscopic streaming term at 
the location of molecule i, u(r^), [|7|], 



m,- 



(9) 



The present treatment uses Pj in the lab frame. A discussion 
of translating CV and its relationship to the peculiar momen- 
tum is given in AppendixlAl 

Finally, the energy density at a point in space is defined by 



N 



p(r,t)S(r,t) = J2{^( r i- r )>f 



i=l 



(10) 



■th 



where the energy of the i molecule is defined as the sum of 
the kinetic energy and the inter-molecular interaction poten- 
tial 6n, 



surface 



U dS. 



(3) 



The rate of change of energy in a CV is expressed in terms of 
fluxes, the pressure tensor and a heat flux vector q, 



lJ v p£dV=-f s [ P S U . 



U-u + q]-dS, 



(4) 



here the energy change due to body forces is not included. 
The divergence theorem relates surface fluxes to the diver- 
gence within the volume, for a variable A, 



A-dS 



V-AdV 



(5) 



V 



In addition, the differential form of the flow equations can be 
recovered in the limit of an infinitesimal control volume 113511 . 



A = lim — 

V~>0 V 



A-dS. 



(6) 



B. Relationship Between the Continuum and the Mesoscopic 
Descriptions 

A mesoscopic description is a temporal and spatial average 
of the molecular trajectories, expressed in terms of a proba- 
bility function, / Irving and Kirkwood ||8|] established the 
link between the mesoscopic and continuum descriptions us- 
ing the Dirac 6 function to define the macroscopic density at 
a point r in space, 



iV 



«=i 



/°( r >*) -Yl \ m i 5 ( r i - r );/ 



(7) 



The angled brackets (a; f) denote the inner product of a with 
/, which gives the expectation of a for an ensemble of sys- 
tems. The mass and position of a molecule i are denoted mj 



ev- 



Pi 
2rrii 



1 N 

-Y 



■'./ 



(ID 



It is implicit in this definition that the potential energy of an 
interatomic interaction, <fiij, is divided equally between the 
two interacting molecules, i and j. 

As phase space is bounded, the evolution of a property, a, 
in time is governed by the equation, 



l ( 



N 



da Pj da 



<9pj mi dr { I 



(12) 



where Fj is the force on molecule i, and a = a(rj(i), Pj(i)) 
is an implicit function of time. Using Eq. ([121) . Irving and 
Kirkwood [8] derived the time evolution of the mass (from 
Eq. |7]), momentum density (from Eq. [8]) and energy density 
(from Eq. ITOb for a mesoscopic system. A comparison of the 
resulting equations to the continuum counterpart provided a 
term-by-term equivalence. Both the mesoscopic and contin- 
uum equations were valid at a point; the former expressed in 
terms of Dirac 6 and the latter in differential form. In the 
current work, the mass and momentum densities are recast 
within the CV framework which avoids use of the Dirac S 
functions directly, and attendant problems with their practi- 
cal implementation. 

III. THE CONTROL VOLUME FORMULATION 

In order to cast the governing equations for a discrete 
system in CV form, a 'selection function' i?j is introduced, 
which isolates those molecules within the region of interest. 
This function is obtained by integrating the Dirac 6 func- 
tion, <5(rj — r), over a cuboid in space, centered at r and 
of side length Ar as illustrated in figure [T(a)| J37J1. Using 
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FIG. 1. (Color online) The CV function and its derivative applied to a system of molecules. The figures were generated using the VMD 
visualization package, 1 36]. From left to right, (a) Schematic of $i which selects only the molecules within a cube, (b) Location of cube 
center r and labels for cube surfaces, (c) Schematic of d$i/dx which selects only molecules crossing the x + and x~~ surface planes. 



S(ri-r)-- 
integral is, 



8{xi — x)5(yi — y)5(zi — z), the resulting triple 



■&i - S ( x i ~ x ) S (Vi ~ V) s ( z i - z)dxdydz 

x~ y~ z~ 



H(xi - x)H(y t - y)H(zi - z) 

= [H{x + - Xi ) - H(x~ - a*)] 
x [H(y+ - Vi ) - H(y- - Vi )] 
x [H{z + - zi) - H(z~ - Zi )] , 



(13) 



where H is the Heaviside function, and the limits of integra- 
tion are defined as, r = r — ^f and r + = r + Qf, for each 
direction (see Fig. |l(b)| i. Note that #j can be interpreted as 
a Lagrangian-to-Control-Volume conversion function {CCV) 
for molecule i. It is unity when molecule i is inside the 
cuboid, and equal to zero otherwise, as illustrated in Fig. 
| l(a)| Using L'Hopital's rule and defining, Ay = AxAyAz, 
the CCV function for molecule i reduces to the Dirac S func- 
tion in the limit of zero volume, 



<J(r 



lim TT7- 

AV->0 AV 



The spatial derivative in the x direction of the CCV function 
for molecule i is, 



dx dxj 



Xi) - S(x - Xi)] S x 



(14) 



where S T i is 



S xi = [H(y+ - Vl ) - H(y~ - Vi j\ 
[H(z+ - Zi ) - H(z~ - Zi)] 



(15) 



Eq. ( fT4l > isolates molecules on a 2D rectangular patch in the 
yz plane. The derivative ddi/dx is only non-zero when 
molecule i is crossing the surfaces marked in Fig. |l(c)| nor- 
mal to the x direction. The contribution of the i molecule 
to the net rate of mass flux through the control surface is ex- 
pressed in the form, Pj • <iSj. Defining for the right x surface, 



dS x h i = S(x + ~xi)S xi , 



(16) 



and similarly for the left surface, dS ,, the total flux Eq. (TT4b 
in any direction r is then, 

dfl- 

-^ = dS+-dSr = dSi. (17) 

The CCV function is key to the derivation of a molecular- 
level equivalent of the continuum CV equations, and it will 
be used extensively in the following sections. The approach 
in sections Ull Al IHI Bl and lHI Dl shares some similarities with 
the work of Serrano and Espanol [38] which considers the 
time evolution of Voronoi characteristic functions. However 
the CCV function has precisely defined extents which allows 
the development of conservation equations for a microscopic 
system. In the following treatment, the CV is fixed in space 
(i.e., r is not a function of time). The extension of this treat- 
ment to an advecting CV is made in AppendixlAl 

A. Mass Conservation for a Molecular CV 

In this section, a mesoscopic expression for the mass in a 
cuboidal CV is derived. The time evolution of mass within 
a CV is shown to be equal to the net mass flux of molecules 
across its surfaces. 

The mass inside an arbitrary CV at the molecular scale can 
be expressed in terms of the CCV as follows, 



p(r,t)dV 



N 

E 



= / V/mjJ^ 

Jvfr[\ 

x+y+z + 



r);f)dV 



);f\dxdydz 



x y z 



N 



*=1 



= J2( m ^uf)- (18) 



Taking the time derivative of Eq. (fT8l and using Eq. (fT2l . 



d_ 



J v P(v,t)dV = ^Yl( m ^f 



N 



-E(^-^+ F -|-^/)- (19) 

The term drriidi/dpi = 0, as i?j is not a function of p,j. 
Therefore, 






(20) 
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where the equality, 9i?j/9tj = —d$i/dr has been used. 
From the continuum mass conservation given in Eq. ([TJ, the 
macroscopic and mesoscopic fluxes over the surfaces can be 
equated, 

6 N , . 

E / pu-dS f = J2(Pi-^i-J)- (21) 

faces JS f t=l ^ ' 

The mesoscopic equation for evolution of mass in a control 
volume is given by, 



where the continuum expression {puu x } + is the average flux 
through a flat region in space with area AA^ = AyAz. This 
kinetic component of the pressure tensor is discussed further 
in Section HITCl 

The configurational term of Eq. (l24t is, 



A' 



d 



cr=E\ F ^^TP^//=E( F ^;/ 



A 



(27) 







A 



dt 



i=l 



A 



sEK^ =-E ft"S«:/ 



i=i 



(22) 



where the total force Fj on particle % is the sum of pairwise- 
additive interactions with potential </>jj , and from an external 
potential ip{. 



AppendixlBlshows that the surface mass flux yields the Irving 
and Kirkwood [8[] expression for divergence as the CV tends 
to a point (i.e. V — >■ 0), in analogy to Eq. ©. 

B. Momentum Balance for a Molecular CV 

In this section, a mesoscopic expression for time evolution 
of momentum within a CV is derived. The starting point is to 
integrate the momentum at a point, given in Eq. ©, over the 
CV, 



A' 



Wi 



d_ 
"On 



A 



E &J + & 



J¥* 



It is commonly assumed that the potential energy of an inter- 
atomic interaction, </>jj, can be divided equally between the 
two interacting molecules, i and j, such that, 



N r 



N ai 1 N 



• J J> 



dvi 



Oth 



(28) 



f p(r, t)u(r, t)dV = V (pA;f)- (23) 7 here the notatio " X& = EJIi E& 

J V r— j \ / for conciseness. Therefore, the conhgu 



Following a similar procedure to that in section HIIAI the for- 
mula (TT2l is used to obtain the time evolution of the momen- 
tum within the CV, 



|/ / (r )i Mr, i )dy = ||(p^;/ 



has been introduced 
igurational term can be 



N , . N . 

c r= oE (%<W) + E( f ^xt^;/ 



'■./ 



i=l 



(29) 



A 



i=l 



■E<^-5^ +, '-^ i/ 



(24) 



where f tj = -dfaj/dri = d^ji/drj and f iext = 
—difji/dri. The notation, #y =${ — flj, is introduced, which 
is non-zero only when the force acts over the surface of the 
CV, as illustrated in Fig. [2] 



K T 



-r 



where the terms Kj- and Cj- are the kinetic and configura- 
tional components, respectively. The kinetic part is, 






A 







i=l x * * ' i=l 



A 



'"/ 



9r ?: J 



(25) 



where p^p, is the dyadic product. For any surface of the CV, 
here x + , the molecular flux can be equated to the continuum 
convection and pressure on that surface, 



P(x + , y, z, t)u(x + , y, z, t)u x (x + ,y, z, t)dydz 



A 



K+dydz = E ( ^dS^f 



*=1 



where K+ is the kinetic part of the pressure tensor due to 
molecular transgressions across the x + CV surface. The av- 
erage molecular flux across the surface is then, 






v 

^ 









&*$* 







A 



{pUU x }~ 



K~ 



1tE(^^/: 



AA+ . 
x i=i 






FIG. 2. (Color online) A section through the CV to illustrate the 
role of #y in selecting only the i and j interactions that cross the 
bounding surface of the control volume. Due to the limited range of 
interactions, only the forces between the internal (red) molecules i 
and external (blue) molecules j near the surfaces are included. 
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Substituting the kinetic {Kj<) and configurational (C7 1 ) 
terms, from Eqs. (l25T l and j29] l into Eq. d24b , the time evo- 
lution of momentum within the CV at the mesoscopic scale 



is. 



N 



at 



i=l 



N 



sEp*/=-E?« 



i=l 



iJV . \ N / 



*J 



1=1 



(30) 



Equations ( l22b and (l30b describe the evolution of mass and 
momentum respectively within a CV averaged over an en- 
semble of representative molecular systems. As proposed by 
Evans and Morriss H, it is possible to develop microscopic 
evolution equations that do not require ensemble averaging. 
Hence, the equivalents of Eqs. d22b and (f30T > are derived for 
a single trajectory through phase space in section HV Al inte- 
grated in time in section IIVBI and tested numerically using 
molecular dynamics simulation in section ITV CI 

The link between the macroscopic and mesoscopic treat- 
ments is given by equating their respective momentum Eqs. 
© and C£D, 

- j> puu ■ dS + F surface + F body 

N , . 

N . \ N / \ 

+ ?Ew + E w4 ( 31 > 



'..; 



*=1 



As can be seen, each term in the continuum evolution of mo- 
mentum has an equivalent term in the mesoscopic formula- 
tion. 

The continuum momentum Eq. (0 can be expressed in 
terms of the divergence of the pressure tensor, II, in the con- 
trol volume from, 

— / pudV = ~ f [puu + n] • dS + F body (32a) 
r — • [puu + n] dV + F body . (32b) 

In the following subsection, the right hand side of Eq. OTI ) 
is recast first in divergence form as in Eq. (I32bl ). and then in 
terms of surface pressures as in Eq. (132at . 



C. The Pressure Tensor 

The average molecular pressure tensor ascribed to a con- 
trol volume is conveniently expressed in terms of the CCV 
function. This is shown inter alia to lead to a number of 
literature definitions of the local stress tensor. In the first 
part of this section, the techniques of Irving and Kirkwood 
]8[] are used to express the divergence of the stress (as with 
the right hand side of Eq. (I32bl )) in terms of intermolecu- 
lar force. Secondly, the CV pressure tensor is related to the 
Volume Average (VA) formula ( 11161 \VW ) and, by considera- 
tion of the interactions across the surfaces, to the Method Of 
Planes (MOP) HKM]. Finally, the molecular CV Eq. (0 
is written in analogous form to the macroscopic Eq. (I32al l. 

The pressure tensor, II, can be decomposed into a kinetic 
k term, and a configurational stress er. In keeping with the 



engineering literature, the stress and pressure tensors have 
opposite signs, 



n = k — a. 



(33) 



The separation into kinetic and configurational parts is made 
to accommodate the debate concerning the inclusion of ki- 
netic terms in the molecular stress IM I39J 14011 . 

In order to avoid confusion, the stress, er, is herein de- 
fined to be due to the forces only (surface tractions). This, 
combined with the kinetic pressure term k, yields the total 
pressure tensor II first introduced in Eq. (f3]). 

1. Irving Kirkwood Pressure Tensor 

The virial expression for the stress cannot be applied lo- 
cally as it is only valid for a homogeneous system, lfl2ll . The 
Irving and Kirkwood [8J] technique for evaluating the non- 
equilibrium, locally-defined stress resolves this issue, and is 
herein extended to a CV. To obtain the stress, a, the inter- 
molecular force term of Eq. (Bil l is defined to be equal to the 
divergence of stress, 

N 



i,3 

litf v (^ ^ - r ^- 6( v - r )] > f ) dv - (34) 



Irving and Kirkwood ||8j] used a Taylor expansion of the Dirac 
S functions to express the pair force contribution in the form 
of a divergence, 



% [S(*i - r) - S(rj 



— ■i i jY ij O ij 6{r i -r) 



where r« = r; — Tj, and Oij is an operator which acts on the 
Dirac 5 function, 



n h 1 d 



J_ / _d_ 
n! V^dr, 



n-l 



(35) 



Equation (l34b can therefore be rewritten, 

h3 

OijSir, - r);/W (36) 

The Taylor expansion in Dirac S functions is not straightfor- 
ward to evaluate. This operation can be bypassed by integrat- 
ing the position of the molecule i over phase space [110 . or by 
replacing the Dirac S with a similar but finite-valued function 
of compact support 1151 Il8l Il9l I2III . In the current treatment, 
the CCV function, $, is used, which is advantageous because 
it explicitly defines both the extent of the CV and its surface 
fluxes. The pressure tensor can be written in terms of the 
CCV function by exploiting the following identities (see Ap- 
pendix of Ref. Ji), 



1 
Oij5(ri - r) = 5(r- r, + sr^ds, 



(37) 
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Equation d36l can therefore be written as, 



V or 



1 N I d 

h3 



1 



x 6(r- r, ; + sr i:j )ds;f)dV. (38) 



Equation Eq. ( 1381 leads to the VA and MOP definitions of 
the pressure tensor. 

2. VA Pressure Tensor 

definition of the stress tensor of Lutsko IU6I1 and Cormier 
et al. um can be obtained by rewriting Eq. (l38l ) as, 



dr J v 



N 



i-titi* 



1 -J 



1 

x A«J(r-ri + sry)ds;/\dV. (39) 



Equating the expressions inside the divergence on both sides 
of Eq. ( 1391 , 14111 . and assuming the stress is constant within 
an arbitrary local volume, AV, gives an expression for the 

VA stress, 



VA 
er = 



yj7 / E \ f iJ r v S ( r ~ r i + srij)ds;f\dV. 
hi n 



2AV 



(40) 



Swapping the order of integration and evaluating the integral 
of the Dirac S function over AV gives a different form of the 
CCV function, $ s , 

■& 8 = I 5(r - n + sr i:j )dV = 
Jv 

[H(x + - Xi + sxij) - H(x~ - Xi + sxij)] 
x [H(y + - Vi + syij) - H(y~ - y L + syij)] 
x [H(z + -Zi + sz i: j) - H(z~ -Zi + szij)] , (41) 

which is non-zero if a point on the line between the two 
molecules, r.; — stij, is inside the cubic region (c.f. rj with 
■di). Substituting the definition, § s (Eq. ETJ, into Eq. ( 1401 
gives, 



VA 



\\r 7 , \*ij r ijkjlf 



2AV" 



1 -j 



(42) 



where ly is the integral from r^ (s = 0) to rj (s = 1) of the 
■& s function, 



kj = I finds. 



Therefore, l^ is the fraction of interaction length between i 
and j which lies within the CV, as illustrated in Fig. [3] The 
definition of the configurational stress in Eq. (l42l) is the same 
as in the work of Lutsko |16] and Cormier et al. IY% . The 
microscopic divergence theorem given in Appendix|A]can be 




FIG. 3. (Color online) A plot of the interaction length given by the 
integral of the selecting function # s defined in Eq. <41t along the 
line between n and rj. The cases shown are for two molecules 
which are a) both inside the volume (kj — 1) and b) both outside 
the volume with an interaction crossing the volume, where kj is the 
fraction of the total length between i and j inside the volume. The 
line is thin (blue) outside and thicker (red) inside the volume. 

applied to obtain the volume averaged kinetic component of 
the pressure tensor, Kq-, in Eq. d2~5l l. 



VA VA 

{puu} + K, 



N 

E 



P f-dS. 

i=\ x " l ' i=l 



a N 

O v^ / PiP, 



'^-s-?: 



#i\f 



Note that the expression inside the divergence includes both 

VA 
the advection, {puu}, and kinetic components of the pres- 
sure tensor. The VA form IU7ll is obtained by combining the 

above expression with the configurational stress a , 

VA VA VA VA VA 
{puu} + k — a = {puu} + II 

1 N I l N \ 



AV <—* \ m,- 



1,3 



In contrast to the work of Cormier et al. IU7I1 . the advection 
term in the above expression is explicitly identified, in order 
to be compatible with the right hand side of Eq. ( I32bb and 
definition of the pressure tensor, II. 

3. MOP Pressure Tensor 

The stress in the CV can also be related to the tractions 
over each surface. In analogy to prior use of the molecular 
CCV function, $j, to evaluate the flux, the stress CCV func- 
tion, iD s , can be differentiated to give the tractions over each 
surface. These surface tractions are the ones used in the for- 
mal definition of the continuum Cauchy stress tensor. The 
surface traction (i.e., force per unit area) and the kinetic pres- 
sure on a surface combined give the MOP expression for the 
pressure tensor [fl3ll . 
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In the context of the CV, the forces and fluxes on the six 
bounding surfaces are required to obtain the pressure inside 
the CV. It is herein shown that each face takes the form of 
the Han and Lee I14J1 localization of the MOP pressure com- 
ponents. The divergence theorem is used to express the left 
hand side of Eq. ( l38l l in terms of stress across the six faces 
of the cube. The mesoscopic right hand side of Eq. d38l ) can 
also be expressed as surface stresses by starting with the CCV 
function $.», 



6 

faces f 



cr ■ dS 



1 



1 N I 



'■j 



dr 



ds;f 



The procedure for taking the derivative of i) s with respect to 
r and integrating over the volume is given in Appendix [C] 
The result is an expression for the force on the CV rewritten 
as the force over each surface of the CV. For the x + face, for 
example, this is, 

N 






sgn{x + - Xi)] S+- ■;/ 



The combination of the signum functions and the ST- ■ term 
specifies when the point of intersection of the line between i 
and j is located on the x + surface of the cube (see Appendix 
[Ct . Corresponding expressions for the y and z faces are de- 
fined by S~^- : when a = {y, z} respectively. 

The full expression for the MOP pressure tensor, which 
includes the kinetic part given by Eq. (126*1 ). is obtained by 
assuming a uniform pressure over the x + surface, 



/ II • dSt = [k - er] ■ n+AA+ 
■>s+ 

= [K+-T+]AA+=P+A,4+ 



(44) 



where n+ is a unit vector aligned along the x coordinate axis, 
n+ = [+1,0,0]; T+ is the configurational stress (traction) 
and P+ the total pressure tensor acting on a plane. Hence, 

1 N /-- 

Vi?™ 6{xi _ x + ) S+,f 






N 



x hi ' 

(45) 

where the peculiar momentum, p, has been used as in Todd 
et al. 11311. If the x + surface area covers the entire domain 



(SZjh = 1 in Eq. d45ll). the MOP formulation of the pressure 

j — i 
is recovered IH3H . 

The extent of the surface is defined through ST- ■, in Eq. 
(1451 1 which is the localized form of the pressure tensor con- 
sidered by Han and Lee 111411 applied to the six cubic faces. 
For a cube in space, each face has three components of stress, 
which results in 18 independent components over the total 
control surface. The quantity, 

dS aij = - [sgn(r+ - r aj ) - sgn(r+ - r ai )] 5+^. 
~2 [ s 9 n ( r a ~ r aj) ~ sgn(r~ - r ai j\ S*^-, 
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T> 






°°« 






X 













dp 



^ 

K* 







W 



<4 






FIG. 4. (Color online) Representation of those molecules selected 
through dS X ij in Eq. d46b with molecules i on the side of the surface 
inside the CV (red) and molecules j on the outside (blue). The CV 
is the inner square on the figure. 



selects the force contributions across the two opposite faces; 
similar notation to the surface molecular flux, dS ij = dSf- — 

dS^ (c.f. Eq. (fTTI)). is used. The case of the two x planes 
located on opposite sides of the cube is illustrated in Fig. [4] 

Taking all surfaces of the cube into account yields the final 
form, 



" r 1 / N 

J2 J a - dS / = - 2 £ ( fi i Yl dS ™ff 

faces f i,j a=l ' 

= --J2( f ij il -dSij>f 



I -J 



(46) 



The vector n, obtained in AppendixICl is unity in each direc- 
tion. The tensor ^ is defined, for notational convenience, to 
be the outer product of the intermolecular forces with n, 



«y = - { ij A = ~ { ij i l l !] = - 



Jxij Jxij Jxij 
Jyij Jyij Jyij 
Jzij Jzij Jzij 



In this form, the #j,- function for all interactions over the 
cube's surface is expressed as the sum of six selection func- 
tions for each of the six faces, i.e. #jj = — Y^c\=l dS a ij. 



4. Relationship to the continuum 

The forces per unit area, or 'tractions', acting over each 
face of the CV, are used in the definition of the Cauchy stress 
tensor at the continuum level. For the x + surface, the traction 
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vector is the sum of all forces acting over the surface, 



1 N I 



4AAJ 



i -j 



sgn(x + ~ Xl )]S+ j; f\ (47) 



which satisfies the definition, 



r± 



a n 



± 



of the Cauchy traction 142H . A similar relationship can be 
written for both the kinetic and total pressures, 



K x — K ' n x J 



± 



n n 



where n~ is a unit vector, n^F = f±l 



01 



The time evolution of the molecular momentum within a 
CV ( Eq. (130b). can be expressed in a similar form to the 
Navier-Stokes equations of continuum fluid mechanics. Di- 
viding both sides of Eq. (l30b by the volume, the following 
form can be obtained; note that this step requires Eqs. Eq. 
Eq. (g5j and Eq. ( 1471 : 



1 d 
AV dt 



N 



i=l 



J2\p^f) + 



{pil a Up} + - {pu a UpY 






Ar 
N 



Arr- 



ir rf 



AV 



i=\ 



A y 2-/ \faiext®i>f 



where index notation has been used (e.g. T x 



(49) 
T± ) with 



the Einstein summation convention. 

In the limit of zero volume, each expression would be simi- 
lar to a term in the differential continuum equations (although 
the pressure term would be the divergence of a tensor and not 
the gradient of a scalar field as is common in fluid mechan- 
ics). The Cauchy stress tensor, er, is defined in the limit that 
the cube's volume tends to zero, so that T + and T~ are re- 
lated by an infinitesimal difference. This is used in contin- 
uum mechanics to define the unique nine component Cauchy 
stress tensor, dcr/dx = lim^ ;c ^.Q[T + + T _ ]/Ax. This limit 
is shown in AppendixlBlto yield the Irving and Kirkwood |8[] 
stress in terms of the Taylor expansion in Dirac 5 functions. 

Rather than defining the stress at a point, the tractions 
can be compared to their continuum counterparts in a fluid 
mechanics control volume or a solid mechanics Finite Ele- 
ments (FE) method. Computational Fluid Dynamics (CFD) 
is commonly formulated using CV and in discrete simula- 
tions, Finite Volume fl. Surface forces are ideal for coupling 
schemes between MD and CFD. Building on the pioneering 
work of O'Connell and Thompson 12311 . there are many MD 
to CFD coupling schemes - see the review paper by Mo- 
hamed and Mohamad 14311 . More recent developments for 
coupling to fluctuating hydrodynamics are covered in a re- 
view by Delgado-Buscalioni 14411 . A discussion of coupling 
schemes is outside the scope of this work, however finite vol- 
ume algorithms have been used extensively in coupling meth- 
ods ll3Ul32Ll45l - l47ll together with equivalent control volumes 
defined in the molecular region. An advantage of the herein 
proposed molecular CV approach is that it ensures conser- 
vation laws are satisfied when exchanging fluxes over cell 



surfaces — an important requirement for accurate unsteady 
coupled simulations as outlined in the finite volume coupling 
of Delgado-Buscalioni and Coveney 14511 . For solid coupling 
schemes, I130ll . the principle of virtual work can be used with 
tractions on the element corners (the MD CV) to give the 
state of stress in the element 04811 . 



a ■ VN a dV = i, N a TdS, 

\ Js 



(50) 



where iV a is a linear shape function which allows stress to 
be defined as a continuous function of position. It will be 
demonstrated numerically in the next section. HVl that the CV 
formulation is exactly conservative: the surface tractions and 
fluxes entirely define the stress within the volume. The trac- 
tions and stress in Eq. (1501 are connected by the weak formu- 
lation and the form of the stress tensor results from the choice 
of shape function N a . 



D. Energy Balance for a Molecular CV 

In this section, a mesoscopic expression for time evolution 
of energy within a CV is derived. As for mass and momen- 
tum, the starting point is to integrate the energy at a point, 
given in Eq. ( 110) . over the CV, 



j P (r,t)£(r,t)dV = Y,(e i # i ;f 



(51) 



The time evolution within the C V is given using formula (112b . 



d_ 



J p{r, t)S(r, t)dV = §- t Y, \ e ^f 



N 



i=l 



f-T \ mi <9r, ; dp,- 



(52) 



Evaluating the derivatives of the energy and CCV function 
results in, 



« N 1 \ 1 N 



dt 



i=l 



i.J 
N 



Pi.f +Pi. f 

rrij rrij 



Vilf 



i=l 



*■ — * \ m ■ m ■ 



111; 



Using the definition of Fj, Newton's 3rd law and relabelling 
indices, the intermolecular force terms can be expressed in 
terms of the interactions over the CV surface, $,■ 



''.;• 



N 



dt 



i=l 



N 



££ («*,;/ =-£ («£■*** 



i.j 



N 



N , v N , , 



«=i 
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The right hand side of this equation is equated to the right 
hand side of the continuum energy Eq. [4] 



energy flux 



heat flux pressure heating 



p£ u ■ dS — f q ■ dS — f II • u ■ dS 

s Js Js 



N 



8=1 



J2\ e izr- dS i'>f 



1 N / ' 



l-:l 



(53) 



where the energy due to the external (body) forces is ne- 
glected. The tij'&ij has been re-expressed in terms of surface 
tractions, <aj ■ dSij, using the analysis of the previous sec- 
tion. In its current form, the microscopic equation does not 
delineate the contribution due to energy flux, heat flux and 
pressure heating. To achieve this division, the notion of the 
peculiar momentum at the molecular location, u(rj) is used 
together with the velocity at the CV surfaces u(r^), follow- 
ing a similar process to Evans and Morriss H. 

IV. IMPLEMENTATION 

In this section, the CV equation for mass, momentum and 
energy balance, Eqs. (l22b . (f30b and ( T53l . will be proved to ap- 
ply and demonstrated numerically for a microscopic system 
undergoing a single trajectory through phase space. 

A. The Microscopic System 

Consider a single trajectory of a set of molecules through 
phase space, defined in terms of their time dependent coor- 
dinates Tj and momentum p i . The CCV function depends on 
molecular coordinates, the location of the center of the cube, 
r, and its side length, Ar, i.e., #j = $j(rj(t),r, Ar), The 
time evolution of the mass within the molecular control vol- 
ume is given by, 



d 

Jt 



N N 

N 
i=l 



i=l 

N 






i=l 



dt Ota 



dt 

dSi, 



(54) 



using, Pj = midri/dt. The time evolution of momentum in 
the molecular control volume is, 



d 



N 



^X>^fo(*)> r ' Ar ) 



dt 



N 



= E 

«=i 

As, dpjdt = Fj, then 

A? 



i=l 

N 

E 

i=l 

dr. 



Il>*. 



dt 



N 



1 dt 

dr; 



dt 



dPi 



dt 



dt 



i=l 



i=l 



PjPi 

mi 



N 

-E 
i=i 



rrii 



■ dSi 



E 

:=i 



A 



dSi + Frfi 

N 



,ih 



(55) 



i.j 



i=i 



where the total force on molecule i has been decomposed into 
surface and 'external' or body terms. The time evolution of 
energy in a molecular control volume is obtained by evaluat- 
ing, 



d 



N 



A 



s !>•>< = £ 



dt 



4=1 



1=1 



d'&i de- L 
i ~dt + ~dt 



N N . 

., mi f— ' mi 

i=i i=i 



1 



N r 

E 

'■j 



i ji 



using, dpi/dt = Fj and the decomposition of forces. The 
manipulation proceeds as in the mesoscopic system to yield, 



a N 

IE"* 



A' 



dt 



»=i 



t=i 

N 



Yen^-dS 



N N 

i Y £- ■ i u-di ■ + Y — ■ i 

2 ^— ' rrii ^— ' rrii 

1,3 i=l 



?9- 



(56) 



The average of many such trajectories defined through Eqs. 
( 1541 ), ( 1331 ) and (T56b gives the mesoscopic expressions in Eqs. 
(1221 . d30l and (T53l , respectively. In the next subsection, the 
time integral of the single trajectory is considered. 

B. Time integration of the microscopic C V equations 

Integration of Eqs. d54"l i, (1331 1 and (TSoT l over the time inter- 
val [0, t] enables these equations to be usable in a molecular 
simulation. For the conservation of mass term, 



N T . N 

Y m [Mr) - MO)} = -/!>*■ dSidt (57) 

1=1 Q 1 = 1 

The surface crossing term, dSj, defined in Eq. (TToT l, involves 
a Dirac 6 function and therefore cannot be evaluated directly. 
Over the time interval [0, r], molecule i passes throu gh a 
given x position at times, t X i f~, where k = 1, 2, ..., Nt x 114911 
. The positional Dirac S can be expressed as, 



N t x 



6(xi(t)-x) = Y 



k=l 



5(t - t x i^k) 

\ x i\^xi,k)\ 



(58) 



where \ii{t x i j.)| is the magnitude of the velocity in the x 
direction at time t x i ^. Equation Eq. (158b is used to rewrite 
dSi in Eq. ( 1571 ) in the form, 



dS, 



ai,k : 



S 9 n i^k- T ) - s 9<^i,k -°) 



Sni.k.V'm.k) 



S ai,k( t c 



(59) 



where a = {x, y, z}, and the fluxes are evaluated at times, 
t ■ , and t~- , for the right and left surfaces of the cube, re- 

oa,k ca,k fc ' 

spectively. Using the above expression, the time integral in 
Eq. (1571 ) can be expressed as the sum of all molecule cross- 
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ings, Nt — Nt x + Nt y + Nt z over the cube's faces, 

Accumulation 
* , 

A N t 



A 



Po 



X m i t^( T ) - ^(°)1 = ~ X) X] mi X] T. l± l dS ai,k 



i=l 



i=l fc=l a- 



1\ l^° 



Advection 



(60) 



In other words, the mass in a CV at time t = t minus its 
initial value at t = is the sum of all molecules that cross its 
surfaces during the time interval. 

The momentum balance equation Eq. d55l l. can also be 
written in time-integrated form, 

AT 



5>i(r)0i(r)-Pi(O)i?i 



i=l 



T 


"a A A 




-I 


Ef-^-sEMti-EM 


dt, 





i=l ij i=l 




and using identity (|59l). 


Accumulation Advection 


N A N t 3 


X lp<W*iW-Pi(o)^(o)] +X X p*X r^i dS ai,k 

i=l i=lk=l a=l IPQi| 


TV ^ AT ^ 




= £ / TijWijMdt + X / «i«t(*)^(*)d« ■ 

ij i=1 



Forcing 



(61) 



The integral of the forcing term can be rewritten as the sum, 

/ hj{t)$ij{t)dt ~AtJ2 f ij (*n) #ij (*n) , 



n=l 



where 7Y T is the number time steps. Equation doTI ) can be 
rearranged as follows, 



A 

X 



Pa«(T)l?iC-)-Pa*(O)0i(O) 



i=l 
{pu a up} + - {pwo-u^}" 



-Ay 



Ar 



T a(3 ~ T a(: 



K af5 ~ K af: 

Ar p 



A A T 



A 



r P 



N T AV 



XX^ext^")^^)' ( 62 ) 



i=l n=l 



where the overbar denotes the time average. The time- 
averaged traction in (l62b is given by, 



J a/9 



A A- 



11 

^Ia77 X X faijit^dS^Atn), 



N T 4AA 



i,j ra=l 



The time-averaged kinetic surface pressure in (1621 ) is, 

A N t 



TT± 1 1 



K 



a/3 



Pai(tk)Pj3i(tk) 



2AA $ ^ i \p^(t k )\ ^ kik) 

-{pu Q u /3 } ± . 



The Eq. (1621 demonstrates that the time average of the fluxes, 
stresses and body forces on a CV during the interval to r, 
completely determines the change in momentum within the 
CV for a single trajectory of the system through phase space 
(i.e. an MD simulation). The time evolution of the micro- 
scopic system, Eq. ( l62l . can also be obtained directly by 
evaluating the derivatives of the mesoscopic expression (|49| i 
and invoking the ergodic hypothesis, hence replacing (a;f) 
with i JZ adt. The use of the ergodic hypothesis is justified 
provided that the time interval, r, is sufficient to ensure phase 
space is adequately sampled. 

Finally, there are no new techniques required to integrate 
the energy Eq. l56l 



A 



J2 [ei(r)Mr) - ^(0)^(0)] 



i=X 







A 



A 



M-*-E£-v. 



dt (63) 



i=l z,j 

which gives the final form, written without external forcing, 

Accumulation Advection 
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N N t 3 



Y / [e l (r)Mr)-e t (0)M0)}+J2Y, e ^T^l dS ^ 
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■ 

M*)0ij(*)d«- 



ly /"Pi(*) 



m o 



Forcing 



(64) 



As in the momentum balance equation, the integral of the 
forcing term can be approximated by the sum, 



PiW 



tij(t)dij{t)dt 



AtJ2^^--{ ij (t n )^ ij (t n ), 



n=l 



where N T is the number time steps. 

In the next section, the elements, Accumulation, Advec- 
tion and Forcing in the above equations are computed indi- 
vidually in an MD simulation to confirm Eqs. ( l60t , (loTT l and 
( l64b numerically. 

C. Results and Discussion 

Molecular Dynamics (MD) simulations in 3D are used in 
this section to validate numerically, and explore the statisti- 
cal convergence of, the CV formalism for three test cases. 
The first investigation was to confirm numerically the con- 
servation properties of an arbitrary control volume. The sec- 
ond simulation compares the value of the scalar pressure ob- 
tained from the molecular CV formulation with that of the 
virial expression for an equilibrium system in a periodic do- 
main. The final test is a Non Equilibrium Molecular Dy- 
namics (NEMD) simulation of the start-up of Couette flow 
initiated by translating the top wall in a slit channel geom- 
etry. The NEMD system is analyzed using the CV expres- 
sions Eqs. (l60l , doTI ) and d64b , and the shear pressure was 
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computed by the VA and CV routes. Newton's equations of 
motion were integrated using the half-step leap-frog Verlet 
algorithm, IpOll . The repulsive Lennard-Jones (LJ) or Weeks- 
Chandler- Anderson (WCA) potential lIBITl . 



#(ri,-)=4e 



ij 



12 



£ 



ij 



e, nj < r c , (65) 



was used for the molecular interactions, which is the 
Lennard-Jones potential shifted upwards by e and truncated 
at the minimum in the potential, rij = r c = 2 1 '®£. The 
potential is zero for ?*« > r c . The energy scale is set by e, 
the length scale by I and molecular mass by m. The results 
reported here are given in terms of £, e and m. A timestep of 
0.005 was used for all simulations. The domain size in the 
first two simulations was 13.68, which contained N = 2048 
molecules, the density was p = 0.8 and the reduced tem- 
perature was set to an initial value of T — 1.0. Test cases 
1 and 2 described below are for equilibrium systems, and 
therefore did not require thermostatting. Case 3 is for a non- 
equilibrium system and required removal of generated heat, 
which was achieved by thermostatting the wall atoms only. 

1. Case 1 

In case 1, the periodic domain simulates a constant energy 
ensemble. The separate terms of the integrated mass, mo- 
mentum and energy equations given in d6Qb . (l6Tb and ( l64l ) 
were evaluated numerically for several sizes of CV. The mass 
conservation can readily be shown to be satisfied as it simply 
requires tracking the number of molecules in the CV. The 
momentum and energy balance equations are conveniently 
checked for compliance at all times by evaluating the resid- 
ual quantity, 

Residual = Accumulation — Forcing + Advection, (66) 

which must be equal to zero at all times for the CV equations 
to be satisfied. This was demonstrated to be the case, as may 
be seen in Figs. |5(a)| and |5(b)] for a cubic CV of side length 
1.52 in the absence of body forces. The evolution of momen- 
tum inside the CV is shown numerically to be exactly equal 
to the integral of the surface forces until a molecule crosses 
the CV boundary. Such events give rise to a momentum flux 
contribution which appears as a spike in the Advection and 
Accumulation terms, as is evident in Fig. |5(a)| The residual 
nonetheless remains identically zero (to machine precision) 
at all times. The energy conservation is also displayed in 
Fig. |5(b)| The average error over the period of the simulation 
(100 MD timeunits) was less than 1%, where the average er- 
ror is defined as the ratio of the mean \Residual\ to the mean 
\Accumulation\ over the simulation. The error is attributed 
to the use of the leapfrog integration scheme, a conclusion 
supported by the linear decrease in error as timestep At — > 0. 

2. Case 2 

As in case 1, the same periodic domain is used in case 2 
to simulate a constant energy ensemble. The objective of this 
exercise is to show that the average of the virial formula for 
the scalar pressure, H v i r , applicable to an equilibrium peri- 
odic system, 
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(67) 
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FIG. 5. The various components in Eq. [66] 'Accumulation' ( — ), 
the time integral of the surface force, 'Forcing' ( X ), and momen- 
tum flux term, Advection' ( ) are shown. 'Forcing' symbols are 

shown every 4th timestep for clarity and the insert shows the full 
ordinate scale over the same time interval on the abscissa. From 
top to bottom, (a) Momentum Control Volume, (b) Energy Control 
Volume. 

arises from the intermolecular interactions across the periodic 
boundaries 111211 . The CV formula for the scalar pressure is, 



1 



n cv =-(p+ + 



, p + - 

1 yy 



P, 



.</.</ 



■PJz+Pzz) 



(68) 



where the P^ a normal pressure is defined in Eq. (145b and 
includes both the kinetic and configurational components 
on each surface. Both routes involve the pair forces, fy. 
However, the CV expression which uses MOP counts only 
those pair forces which cross a plane while VA (Virial) sums 
fij r ij over m e whole volume. It is therefore expected that 
there would be differences between the two methods at short 
times, converging at long times. A control volume the same 
size as the periodic box was taken. The time averaged control 
volume, (Jfcv) an d virial (II„j r ) pressure values are shown 
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FIG. 6. Tlvir and Ylcv from Eqs. (J67J and d68t respectively. The 
configurational and kinetic pressures are separated with configura- 
tional values typically having greater magnitudes (~ 4.0) than ki- 
netic (~ 0.6). Continuous lines are control volume pressures and 
dotted lines are virial pressure. 



in Fig.|6]to converge towards the same value with increasing 
time. The simulation is started from an FCC lattice with a 
short range potential (WCA) so the initial configurational 
stress is zero. It is the evolution of the pressure from this 
initial state that is compared in Fig. [6] The virial kinetic 
pressure makes use of the instantaneous values of the domain 
molecule's velocities at every time step. In contrast, the 
CV kinetic part of the pressure is due to molecular surface 
crossings only, which may explain its slower convergence 
to the limiting value than the kinetic part of the virial 
expression. To quantify this difference in convergence for 
the two measures of the pressure, the standard deviation, 
SD(.t), is evaluated, ensuring decorrelation IMTII using block 
averaging pill . For the kinetic virial, SD(n v i r ) = 0.0056, 
and configurational, SD(a v i r ) = 0.0619. For the kinetic 
CV pressure SD(kcv) = 0.4549 and configurational 
SD((Tqy) = 0.2901. The CV pressure, which makes 
use of the MOP formula, would therefore require more 
samples to converge to a steady state value. However, the 
MOP pressures are generally more efficient to calculate 
than the VA. More usefully, from an evaluation of only the 
interactions over the outer CV surface, the pressure in a 
volume of arbitrary size can be determined. 
Figure|7]is a log-log plot of the Percentage Discrepancy (PD) 
between the two (PD = [100 x \U CV - Tl vir \/Tl vir ]). 
After 10 million timesteps or a reduced time of 5 x 10 4 , 
the percentage discrepancy in the configurational part has 
decreased to 0.01%, and the kinetic part of the pressure 
matches the virial (and kinetic theory) to within 0.1%. The 
total pressure value agrees to within 0.1% at the end of this 
averaging period. The simulation average temperature was 
0.65, and the kinetic part of the CV pressure was statis- 
tically the same as the kinetic theory formula prediction, 
K CV — P^bT = 0.52 15111 . The VA formula for the pressure 
in a volume the size of the domain is by definition formally 
the same as that of the virial pressure. The next test case 
compares the CV and VA formulas for the shear stress in a 



FIG. 7. The percentage relative difference between the virial and 
control volume time-accumulated scalar pressures (PD defined in 
the text). Values for the kinetic, configurational and total PD are 
shown. 

system out of equilibrium. 

3. Case 3 

In this simulation study, Couette flow was simulated by 
entraining a model liquid between two solid walls. The top 
wall was set in translational motion parallel to the bottom 
(stationary) wall and the evolution of the velocity profile to- 
wards the steady-state Couette flow limit was followed. The 
velocity profile, and the derived CV and VA shear stresses are 
compared with the analytical solution of the unsteady diffu- 
sion equation. Four layers of tethered molecules were used 
to represent each wall, with the top wall given a sliding ve- 
locity of, Uq = 1.0 at the start of the simulation, time t = 0. 
The temperature of both walls was controlled by ap plyi ng the 
Nose-Hoover (NH) thermostat to the wall atoms [52]. The 
two walls were thermostatted separately, and the equations 
of motion of the wall atoms were, 



mi 

iext= r i ( 4fc4r ?) +6fc6r i)' 
1 



Qn 



■ N _ _ 

y h^p^ - 3T 



(69a) 

(69b) 
(69c) 

(69d) 



where nj is a unit vector in the x— direction, m n = m, and 
f j . is the tethered atom force, using the formula of Petravic 
and Harrowell JH (fc 4 = 5 x 10 3 and (: 6 = 5x 10 6 ). The 
vector, rj„ = r, — rrj, is the displacement of the tethered atom, 
i, from its lattice site coordinate, ro. The Nose-Hoover ther- 
mostat dynamical variable is denoted by £, To = 1.0 is the 
target temperature of the wall, and the effective time constant 



or damping coefficient, in Eq. ( 169dl i was given the value, 
Q^ = NAt. The simulation was carried out for a cubic do- 
main of sidelength 27.40, of which the fluid region extent 
was 20.52 in the y— direction. Periodic boundaries were used 
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FIG. 8. (Color online) Schematic diagram of the NEMD simulation 
geometry consisting of a sliding top wall and stationary bottom wall, 
both composed of tethered atoms. The simulation domain contained 
a lattice of contiguous C V used for pressure averaging (shown by the 
small boxes) while the thicker line denotes a single CV containing 
the entire liquid region. 

in the streamwise (x) and spanwise (z) directions. The re- 
sults presented are the average of eight simulation trajecto- 
ries starting with a different set of initial atom velocities. The 
lattice contained 16384 molecules and was at a density of 
p = 0.8. The molecular simulation domain was sub-divided 
into 4096 (16 3 ) control volumes, and the average velocity 
and shear stress was determined in each of them. A larger 
single CV encompassing all of the liquid region of the do- 
main, shown bounded by the thick line in Fig. [8] was also 
considered. 

The continuum solution for this configuration is consid- 
ered now. Between two plates, there are no body forces and 
the flow eventually becomes fully developed, 15411 so that Eq. 
(0 can be simplified and after applying the divergence theo- 
rem from Eq. (0 it becomes, 

|/^ = -jVrw, 

which is valid for any arbitrary volume in the domain and 
must be valid at any point for a continuum. The shear pres- 
sure in the fluid, H xy (y), drives the time evolution, 



dpu x 
dt 



dU 



J-H 



dy 



For a Newtonian liquid with viscosity, p, |5 

du x 



n. 



-/'■- 



this gives the ID diffusion equation, 

du x p d 2 u x 



dt 



p dy 2 



(70) 



(71) 



assuming the liquid to be incompressible. This can be solved 
for the boundary conditions, 



u x (0,t)=0 



,{L,t) = U u x (y,0) = 0, 
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FIG. 9. The y— dependence of the streaming velocity profile at 
times t = 2 n for n = 0, 2, 3, 4, 5, 6 from right to left. The squares 
are the NEMD CV data values and the analytical solution to the 
continuum equations of Eq. J72t is given at the same six times as 
continuous curves. 

where the bottom and top wall-liquid boundaries are at y = 
and y = L, respectively. The Fourier series solution of these 
equations with inhomogeneous boundary conditions 115511 is, 

U y = L 

MV, t) = \ Yl u n(t)sin (^) < y < L (72) 

y = 



n=l 



where X n = (nn/L) 2 and u n {t) is given by, 



u n {t) 



2E/ (-l) n 



exp 



\nPt 
P 



The velocity profile resolved at the control volume level 
is compared with the continuum solution in Fig. [9] There 
were 16 cubic NEMD CV of side length 1.72 spanning the 
system in the y direction, with each data point on the figure 
being derived from a local time average of 0.5 time units. 
The analytic continuum solution was evaluated numerically 
from Eq. (l72l with n = 1000 and p = 1.6, the latter a 
literature value for the WCA fluid shear viscosity at p = 0.8 
and T = 1.0, Il56ll . There is mostly very good agreement 
between the analytic and NEMD velocity profiles at all 
times, although some effect of the stacking of molecules 
near the two walls can be seen in a slight blunting of the 
fluid velocity profile very close to the tethered walls (located 
by the horizontal two squares on the far left and right of the 
figure) which is an aspect of the molecular system that the 
continuum treatment is not capable of reproducing. 

The VA and CV shear pressure, given by Eqs. (l43l l 
and d45l l. are compared at time t = 10 in Fig. [10] The 
comparison is for a single simulation trajectory resolved 
into 16 cubic volumes of size 1.72 in the y— direction, with 
averaging in the x and z directions and over 0.5 in reduced 
time. The figure shows the shear pressure on the faces of the 
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FIG. 10. The y— dependence of the shear pressure at t = 10, aver- 
aged over 100 timesteps and for a single simulation trajectory. The 
VA value from Eq. d43t are the squares. The CV surface traction 
from Eq. ( 145 1 is indicated by x and o for the top and bottom sur- 
faces, respectively. The solid gray line displays the resulting pres- 
sure field using Eq. d50l l with linear shape functions. 



FIG. 1 1. As Fig. [10] except that the NEMD results are averaged over 
a set of eight independent simulations of 1, 000 timesteps (5 reduced 
time units) each. The simulation-derived VA and CV shear pres- 
sures are compared with the continuum analytical solution given in 
Eq. J73b (solid black line). The jump in the profile on the right of 
the figure is due to the presence of the tethered wall. 



CV. Inside the CV, the pressure was assumed to vary linearly, 
and the value at the midpoint is shown to be comparable 
to the VA-determined value. Figure [10] shows that there 
is good agreement between the VA and CV approaches. 
Note that the CV pressure is effectively the MOP formula 
applied to the faces of the cube, and hence this case study 
demonstrates a consistency between MOP and VA. We have 
shown previously that this is true for the special case of an 
infinitely thin bin or the limit of the pressure at a plane 02211 . 
Practically, the extent of agreement in this exercise is limited 
by the inherent assumptions and spatial resolution of the two 
methods; a single average over a volume is required for VA, 
but a linear pressure relationship is assumed for CV to obtain 
the pressure tensor value corresponding to the center of the 
CV. 

The continuum analytical xy pressure tensor component 
can be derived analytically using the same Fourier series ap- 
proach for du x /dy,$5M, 



Kxy(y,t) = - 



L 



i + 2£(-iy 



p cos (rr-) 



which is valid for the entire domain < y < L. 

A statistically meaningful comparison between the CV, 
VA and continuum analytic shear pressure profiles requires 
more averaging of the simulation data than for the streaming 
velocity, 15711 . and eight independent simulation trajectories 
over 5 reduced time units were used. Figure [TTJ shows 
that the three methods exhibit good agreement within the 
simulation statistical uncertainty. 

As a final demonstration of the use of the CV equations, 
the control volume is now chosen to encompass the entire 
liquid domain (see Fig. [8]), and therefore the external forces 



arise from interactions with the wall atoms only. The mo- 
mentum equation, Eq. ( 1331 ), is written as, 



N 



N 



Q> 



Of 



i=l 

N 



i=l 



Ep^ = -E^-^+E f w^ 



i=l 



y yiijUdxij + Iy UOyij -\~ lij(Xb Z ij J , 



l-J 



which can be simplified as follows. For term, in the 
above equation, the fluxes across the CV boundaries in the 
streamwise and spanwise directions cancel due to the peri- 
odic boundary conditions. Fluxes across the xz boundary 
surface are zero as the tethered wall atoms prevent such cross- 
ings. The force term, (2), also vanishes because across the 
periodic boundary, fijdS^,- = —^ijdS~:-, (similarly for z). 
The external force term, (3), is zero because all the forces 
in the system result from interatomic interactions. The sum 
of the fyij force components across the horizontal bound- 
aries will be equal and opposite, and by symmetry the two 
(73) fzij terms in @ will be zero on average. The above equation 



IZl] 

therefore reduces to, 



d N 



i N 



xi] aJ yij 



Jxijdb, 



!)i J 



(74) 



'■•J 



As the simulation approaches steady state, the rate of change 
of momentum in the control volume tends to zero because 
the difference between the shear stresses acting across the 
top and bottom walls vanishes. The forces on the xz plane 
boundary and momentum inside the CV are plotted in Fig. 
[T2"lto confirm Eq. (1741 numerically. The time evolution of 
these molecular momenta and surface stresses are compared 
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FIG. 12. The evolution of surface forces and momentum change 
for a molecular CV from Eq. 1 1741 . (points) and analytical solution 
for the continuum (Eqs. ((77}, d78t and d76t). presented as lines on 
the figure. The Residual, defined in Eq. d66t . is also given. Each 
point represents the average over an ensemble of eight independent 
systems and 40 timesteps. 

to the analytical continuum solution for the CV, 



d-tJ v PU * dV 



ST J JSJ J 



(75) 



The normal components of the pressure tensor are non-zero 
in the continuum, but exactly balance across opposite CV 
faces, i.e. ELJ^ — ^xx- By appropriate choice of the gauge 
pressure, H xx does not appear in the governing Eq. (l75l l. The 
left hand side of the above equation is evaluated from the an- 
alytic expression for u x , 



— / pu x dV 
dt J v 



2AzA*^ £[!-(-!)' 



P . 



(76) 



The right hand side is obtained from the analytic continuum 
expression for the shear stress, for the bottom surface at y = 
0, 



boundaries. The average response nevertheless agrees well 
with the analytic solution, bearing in mind the element of 
uncertainty in the matching state parameter values. This ex- 
ample demonstrates the potential of the C V approach applied 
on the molecular scale, as it can be seen that computation of 
the forces across the CV boundaries determines completely 
the average molecular microhydrodynamic response of the 
system contained in the CV. In fact, the force on only one of 
the surfaces is all that was required, as the force terms for 
the opposite surface could have been obtained from Eq. (|74| ). 



V. CONCLUSIONS 

In analogy to continuum fluid mechanics, the evolution 
equations for a molecular systems has been expressed within 
a Control Volume (CV) in terms of fluxes and stresses across 
the surfaces. A key ingredient is the definition and manipula- 
tion of a Lagrangian to Control Volume conversion function, 
i9, which identifies molecules within the CV. The final ap- 
pearance of the equations has the same form as Reynolds' 
Transport Theorem applied to a discrete system. The equa- 
tions presented follow directly from Newton's equation of 
motion for a system of discrete particles, requiring no ad- 
ditional assumptions and therefore sharing the same range of 
validity. 

Using the CCV function, the relationship between Volume 
Average (VA) d [H and Method Of Planes (MOP) pres- 
sure II 13l Q4I1 has been established, without Fourier transfor- 
mation. The two definitions of pressure are shown numer- 
ically to give equivalent results away from equilibrium and, 
for homogeneous systems, shown to equal the virial pressure. 

A Navier-Stokes-like equation was derived for the evo- 
lution of momentum within the control volume, expressed 
in terms of surface fluxes and stresses. This pro- 
vides an exact mathematical relationship between molecular 
fluxes/pressures and the evolution of momentum and energy 
in a CV. Numerical evaluations of the terms in the conserva- 
tion of mass, momentum and energy equations demonstrated 
consistency with theoretical predictions. 

The CV formulation is general, and can be applied to de- 
rive conservation equations for any fluid dynamical property 
localised to a region in space. It can also facilitate the deriva- 
tion of conservative numerical schemes for MD, and the eval- 
uation of the accuracy of numerical schemes. Finally, it al- 
lows for accurate evaluation of macroscopic flow properties, 
in a manner consistent with the continuum conservation laws. 



fll xy dS+ = -2AxAz^f^ 

J S f 71=1 
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(77) 



(78) 



In Fig Q~2j the momentum evolution on the left hand side of 
Eq. ( l74l is compared to Eq. d76l) , Equations ( 1771 ) and ( fTHI l are 
also given for the shear stresses acting across the top and bot- 
tom of the molecular control volume (right hand side of Eq. 
(l74l). The scatter seen in the MD data reflects the thermal 
fluctuations in the forces and molecular crossings of the CV 



Appendix A: Discrete form of Reynolds' Transport Theorem 
and the Divergence Theorem 

In this appendix, both Reynolds' Transport Theorem and 
the Divergence Theorem for a discrete system are derived. 
The relationship between an advecting and fixed control vol- 
ume is shown using the concept of peculiar momentum. 

The microscopic form of the continuous Reynolds' Trans- 
port Theorem [Qj] is derived for a property x — x( r jiPi)£) 
which could be mass, momentum or the pressure tensor. The 
CCV function, i?j, is dependent on the molecule's coordinate; 
the location of the cube center, r, and side length, Ar, which 
are all a function of time. The time evolution of the CV is 
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The velocity of the moving volume is defined as u = dr/dt, 
which can be different to the macroscopic velocity u. Sur- 
face translation or deformation of the cube, d'&i/dAr, can be 
included in the expression for velocity u. The above analysis 
is for a microscopic system, although a similar process for a 
mesoscopic system can be applied and includes terms for CV 
movement in Eq. (11211 . 

Hence Reynolds treatment of a continuous medium j^\ is 
extended here to a discrete molecular system, 

, N 



-J2x(m(ri(t),r(t),Ar(t)) 



dt 



N r 

E 
t=i 



dt % 



Pi 

mi 



dS; 



(Al) 



The conservation equation for the mass, % = toj, in a moving 
reference frame is, 



A ; 



N r 



dt 



E m ^ = E 



TO; \ U — 



■ dS.; 



(A2) 



In a Lagrangian reference frame, the translational velocity of 
CV surface must be equal to the molecular streaming veloc- 
ity, i.e., u(r^) = u(rj), so that, 



E 
j=i 



u 



^L ) . dS, 



N 



= - E p * ■ dSi 



The evolution of the peculiar momentum, x — Pj. m a mov- 
ing reference frame is, 



, N N 



dt 
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Here an inertial reference frame has been assumed so that 
dpi/dt = dpi/dt = Fj. For a simple case (e.g. one dimen- 
sional flow) it is possible to utilize a Lagrangian description 
by ensuring, u(r^) = tt(rj), throughout the time evolution. 
In more complicated cases, this is not always possible and 
the Eulerian description is generally adopted. 

Next, a microscopic analogue to the macroscopic diver- 
gence theorem is derived for the generalized function, x> 



N a 

y- 



x( r i,Pi,t)$( r i-r) 



dV 



Jv i=l 



The vector derivative of the Dirac 8 followed by the integral 
over volume results in, 
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where the limits of the cuboidal volume are, r^ 

Ar 



Ar 



and 



r = r — =^ . The mesoscopic equivalent of the continuum 
divergence theorem (Eq. (0) is therefore, 
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Appendix B: Relation between Control Volume and 
Description at a Point 

This Appendix proves that the Irving and Kirkwood |@] 
expression for the flux at a point is the zero volume limit of 
the CV formulation. As in the continuum, the control volume 
equations at a point are obtained using the gradient operator 
in Eq. ([6}, the flux at a point can be shown by taking the zero 
volume limit of the gradient operator of Eq. ©. Assuming 
the three side lengths of the control volume, Ax, Ay and Az, 
tend to zero and hence the volume, AV, tends to zero, 

V • pu = lim lim lim — — - — - — 

Ai->0Ai;->0Az->0 AxAyAz 
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v- / d#i ddi d$ % 

x 2^ \P"'~ + Piy~srr +Piz-^;f 



8=1 



dx 



dy 



dz 



(Bl) 



from Eq. (t2TT> . For illustration, consider the x component 
above, where 



t face 



d-d- ' 

— - = [8(x + - x,j) - S(x~ 



Xi)\S xi . (B2) 



Using the definition of the Dirac 8 function as the limit of two 
slightly displaced Heaviside functions, 



S(n = lim 
A£^0 

the limit of the S T i term is, 



A? 



lim lim S xi = 8{yi - y)8{zi - z) 
Ay— s-0 Az->0 

The Ax —t limit for Xf ace (defined in Eq. (1B2I >) can be 
evaluated using L'Hopital's rule, combined with the property 
of the 8 function, 



d(A0 \} 2 



1 d A it ^ 



so that, 



lim x face =—5(x-Xi). 
Ax->0 J OX 
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Therefore, the limit of ddi/dx as the volume approaches 
zero is, 

T 1- T d^i ^ , . 

hm lim hm — — = — — o (Yi — r) , 
Ai->0Ay->0Az->0 OX OX 

Taking the limits for the x, y and z terms in Eq. (IB lb yields 
the expected Irving and Kirkwood \g\ definition of the diver- 
gence at a point, 
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This zero volume limit of the CV surface fluxes shows that 
the divergence of a Dirac S function represents the flow of 
molecules over a point in space. The advection and kinetic 
pressure at a point is, from Eq. 



y.[ P uu + K } = J^(^-.^S(r i -r);f 



4 = 1 



dr m,-. 



The same limit of zero volume for the surface tractions de- 
fines the Cauchy stress. Using Eq. © and taking the limit of 
Eq. d46b , written in terms of tractions, 
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For the r+ surface, and taking the limits of Ar y and Ar z 
using L'Hopital's rule, 
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The indices /3, m and 7 can be x, y or 2 and f denotes the top 
surface (+ superscript), bottom surface (— superscript) orCV 
center (no superscript). The w selecting function includes 
only the contribution to the stress when the line of interaction 
between i and j passes through the point r' in space. The 
difference between T+ and T~ tends to zero on taking the 
limit Ar x — > 0, so that L'Hopital's rule can be applied. Using 
the property, 
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integral between two molecules introduced in Eq. (1371 ) 
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where the sifting property of the Dirac S function in the r x 
direction has been used to express the integral between two 
molecules in terms of the w xyz function. Hence, 
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As the choice of shifting direction is arbitrary, use of r y or 



in the above treatment would result in zu yzx and m 



J z X y 



re- 



spectively. Therefore, Eq. (1381 1. without the volume integral, 
can be expressed as, 

1 
faijrpij-^— / 5(r-ri + 8Tij)ds;f) 



1 N 

-Y 

h3 



drr- 







1 N 



fi 



'■.;« 



dvo 



X yz 



&W' 



y X z 



dvo 



zxy 



dr x 



dr h 



dr z 



;/ 



As Eq. (I38b is equivalent to the Irving and Kirkwood (8J 
stress of Eq. (136b . the Irving Kirkwood stress is recovered in 
the limit that the CV tends to zero volume. 
This Appendix has proved therefore that in the limit of zero 
control volume, the molecular CV Eqs. (F22b and (|49| i recover 
the description at a point in the same limit that the contin- 
uum CV Eqs. (Q} and (0 tend to the differential continuum 
equations. This demonstrates that the molecular CV equa- 
tions presented here are the molecular scale equivalent of the 
continuum CV equations. 

Appendix C: Relationship between Volume Average and 
Method Of Planes Stress 

This Appendix gives further details of the derivation of the 
Method Of Planes form of stress from the Volume Average 
form. Starting from Eq. (l38b written in terms of the C V func- 
tion for an integrated volume, 
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Taking only the x derivative above, 



then, 
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where G(s) is, 

G(s) = [H(y + - y { + syij) - H(y~ - y t + syij)] 
x [H(z + - Zi + szij) - H{z~~ - Zi + szij)] . 

As 5(ax) = y-rS(x) the XijX~t G(s) term in Eq. (IC21 > can 
be expressed as, 



+ ni \ — 3 



s G(s). (C3) 



The integral can be evaluated using the sifting property of the 
Dirac S function 158fl as follows, 
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where the signum function, sgn(xij) = Xij/\xij\. The S, 
term is the value of s on the cube surface, 
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Repeating the same process for the other faces allows Eq. 
( 1C1I ) to be expressed as, 
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sgn(r ai j)sgn 



[sgn(r± - r aj ) - sgn(r± - r ai )] S^j 
= [111]. This is the force 



over the CV surfaces, Eq. d46*b . in section HIlCI 

To verify the interpretation of SJ-- used in this work, 

XIJ 

consider the vector equation for the point of intersection of 
a line and a plane in space. The equation for a vector a 

between r$ and i\- is defined as a = rj 



— s 



1. 1 



The 



'./ 



plane containing the positive face of a cube is defined by 
(r + p) • n where p is any point on the plane and n is nor- 
mal to that plane. By setting a = p and upon rearrangement 

r, 



of 



.+ 



s y-^-r I ■ n, the value of s at the point of inter- 



section with the plane is, 



The definition ST. . (analogous to S x i in Eq. ( fT3l l) has been 

X% J 

introduced as it filters out those ij terms where the point 
of intersection of line rij and plane x + has y and z com- 
ponents between the limits of the cube surfaces. The cor- 
responding terms, Sfj a , are defined for a — {y, z}. Tak- 
ing H (0) = i, the Heaviside function can be rewritten as 
H(ax) = 2 (sgn(a)sgn(x) — 1), and, 



H 



-H 



= 2 S 9 n [ — ) [sgn(x + - Xj) - sgn(x + - a;,)] , 

so the expression, XijX~t G(s) in Eq. ( 1C2I ) becomes, 

1 
x ij J x jace G ( s ) ds = -^sgn{xij)sgn ( — 



x [sgn(x — Xj) —sgn{x — Xj)]S. 



The signum function, sgn 



XIJ' 



cancels the one obtained 



from integration along s, sgn(xij). The expression for the 
x + face is therefore, 



The point on line a located on the plane is, 
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Taking n as the normal to the x surface, i.e. 
n ^n x . = [1,0,0], then, 
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written using index notation with a = {x, y, z}. The vector 
x+ is the point of intersection of line a with the x + plane. A 
function to check if the point x+ on the plane is located on the 
region between y^ and z , would use Heaviside functions 
and is similar to the form of Eq. (TT31 ), 



^j = [ H (y + -4 P )- H (y~ 

x [H (z + - x+ ) - H {z 



L yp)\ 

4 P )] 



which is the form obtained in the text by direct integration of 
the expression for stress, i.e. Eq. ( IC4I ). 
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